// LennardJones.cpp: archivo de proyecto principal.

#include <complex>
#include <fstream>
#include <iostream>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<time.h>


using namespace std;


#define kB 8.314e-03
#define maxK 10000
#define KsqMax  49
#define Kmax 7
#define pi 3.1415926535898
#define avogadro 602.205

template<typename T>inline T SQR(T x)
	{  return(x*x);  }
template<typename T>inline T POW3(T x)
	{  return(x*x*x);  }

template<typename T>inline T VLENGTH3(const T *x)
	{  return(sqrt(SQR(x[0])+SQR(x[1])+SQR(x[2])));  }


#include "tiempo.h"
#include "utils.h"
#include "particulas.h"
#include "caja.h"
#include "interaccion.h"
#include "termostato.h"
#include "verlet.h"
#include "verletClass.h"
#include "calcValues.h"

//#include "dynamic.h"


int main(void)
{
	//double Temp=0.1;
	double Temp=3.0;
	double delta=1.0;
	double diff=0.005;
	int npart=1000;
	double phi=0.01;
	double rcut=2.5;
	double Etotal;
	double Dt=0.005;

	///Empleo del termostato
	bool termo = true;
	///Reanudar simulacion anterior
	bool restart = false;
	//Impresion de valores energia/presion/temperatura/deriva en pantalla
	bool printScreen = true;
	///Caja de la simulacion
	caja *box;
	///Configuracion para el salvado logaritmico de configuraciones
	int nout=0; ///Salvado de configuraciones
	int nsave=0;///Indice de paso del siguiente salvado de configuraciones
	///Importamos la escala logaritmica de almacenamiento de configuraciones
	int  *NSS;
	NSS=newtempi();
	nsave = NSS[nout];

	char name[50];
	sprintf (name, "%s","H:\\config.xyz");
	ofstream finput (name,ios::out);
	finput.close();
	sprintf (name, "%s","resultados.txt");
	ofstream finp (name,ios::out);
	finp.close();

	sprintf(name,"%s","H:\\config");

	///Comenzamos una simulacion desde el principio?
	if(!restart)
	{
		box = new caja(npart, Temp, phi, rcut,Dt);
	}
	else
	{
		box = new caja("inicio.txt",Dt);
	}

	verlet *vrl= new verlet(box,Dt,termo);

	box->printXYZfile(name, 0);
	//box->listaVecinos();
	box->linkCell();
	interaccionLinkCellList(box);


	//cout<<"nvecinos: "<<box->part[1].nvecinos<<endl;
	for (int i=0;i<1000000;i++)
	{
		//box->followPart(0);
		//if (i%200==0)
			//box->listaVecinos();
		//Etotal=interaccion(box,rcut);
		box->nacc=i;
		//Etotal=interaccionLinkCellList(box);

		//verlet(box,Dt);
		//verletLinkCellList(box,Dt,termo);
		vrl->LeapFrogLinkCellList();
		//vrl->verletLinkCellList();

		if(i==nsave)
		{
			cout<<"Salvando paso "<<i<<endl;
			nout++;
			nsave=NSS[nout];
			printValues(box);
			//cout<<"Etotal: "<<Etotal<<endl;
			//box->printConfig();
			//values(box, i,printScreen);
			box->printXYZfile(name, i);
		}
	}

	delete box;

//	box.Temp=Temp;
//	box.npart = npart;
//	box.initSquareBox(npart, Temp, phi);
//	dynamic din;
//	din.initial(box, delta, diff);
//
//	cout<<"*************** Ewald"<<endl;
//	ewald ewd;
//	ewd.initial(box);
//	ewd.setupEwald();
//	ewd.initialCharges();
//	ewd.ewrForcesCharges();
//	ewd.ewkForcesCharges();
//	ewd.sumForces();
//
//	for(int i=0;i<1000;i++)
//	{
//		cout<<i<<endl;
//	din.box = ewd.box;
//	din.gaussian();
//	din.updateR();
//	din.updateV();
////
////
////	ewd.box = din.box;
////
////	ewd.ewrForcesCharges();
////	ewd.ewkForcesCharges();
////	ewd.sumForces();
//	}


	cout<<"sacabo"<<endl;
	system("PAUSE");
	return EXIT_SUCCESS;
}
